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Abstract. In this paper we examine the issue of characterizing the transport associated with grav- 
itational instabihties in relatively cold discs, discussing in particular under which condition it can 
be described within a local, viscous framework. We present the results of global, three-dimensional, 
SPH simulations of self-gravitating accretion discs, in which the disc is cooled using a simple 
parameterization for the cooling function. Our simulations show that the disc settles in a "self- 
regulated" state, where the axisymmetric stability parameter Q~l and where transport and energy 
dissipation are dominated by self-gravity. We have computed the gravitational stress tensor and com- 
pared our results with expectations based on a local theory of transport. We find that, for disc masses 
smaller than 0.25M* and aspect ratio H/r <0.l, transport is determined locally, thus allowing for a 
viscous treatment of the disc evolution. 



INTRODUCTION 

One of the basic unknowns in accretion disc theory is the physical mechanism ukimately 
responsible for angular momentum transport and energy dissipation in the disc. The 
usual way to overcome this difficulty is to assume that transport is dominated by some 
kind of viscous process, and to give an ad hoc prescription for the component of 
the stress tensor T (usually vertically integrated, under the assumption that the disc is 
geometrically thin). The most widely used prescription is the so-called a -prescription 
yj, according to which: 
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where Q. is the angular velocity of the disc, Cs is the thermal speed, E is the surface 
density, and a a free parameter, expected to be smaller than unity. 

It has recently been recognized that accretion discs threaded by a weak magnetic field 
are subject to MHD instabilities (see Balbus and Hawley [2] and references therein), 
that can induce turbulence in the disc, thereby being able to transport angular momen- 
tum outwards and to promote the accretion process. However, in many astrophysically 
interesting cases, such as the outer regions of protostellar discs, the ionization level is 
expected to be very low, reducing significantly the effects of magnetic fields in deter- 
mining the dynamics of the disc. A possible alternative source of transport in cold discs 
is provided by gravitational instabilities. 



The axisymmetric stability of a thin disc with respect to gravitational disturbances is 
determined by the condition |3]: 

where k is the epicyclic frequency. The external regions of many observed systems 
are likely to be gravitationally unstable. If we consider, for example, the case of AGN 
accretion discs, it can be shown [|4|] that Q falls below unity already at a distance of 
~ 10^ pc from the central black hole. In the case of protostellar discs, models of the 
outburst in FU Orionis systems (a class of young stellar objects that experience a phase 
of enhanced accretion) show that the disc is marginally stable already at a distance of 
fii 1 AU from the central young star 12] . 

It has been suggested [6, 7] that, in a cold enough disc, where condition Q is violated, 
gravitational instabilities would provide a source of energy dissipation, therefore heating 
up the disc and leading to a self-regulated state, with the disc close to marginal stability. 
Gravitationally induced spiral structures are also able to transport angular momentum. 
In this paper we examine the issue of characterizing the angular momentum transport 
and energy dissipation in self-gravitating discs. In particular, we would like to test the 
extent to which the transport due to gravitational instabilities can be described as a local, 
viscous process, amenable to a description analogous to Eq. dl]), as proposed by Lin and 
Pringle 

Balbus and Papaloizou [9] have argued that the long-range nature of the gravitational 
field would preclude a local description of transport, as implicitly assumed in a viscous 
scenario. In particular, they have shown how the energy flux contains extra terms, related 
to wave transport, which are "anomalous" with respect to a viscous flux. Gammie 
ilOl has performed local, shearing sheet, zero-thickness simulations of self-gravitating 
accretion discs, and has concluded that a local description is adequate in such "razor- 
thin" discs, extrapolating his results to discs for which H/r < 0.1. However, Gammie's 
simulations are not appropriate to test global effects, since locality is set up from the 
beginning, and they are only valid for infinitesimally thin discs. Rice et al. [H, 1,J2J have 
already shown using global, 3D simulations, how global effects can be important in the 
dynamics of self-gravitating discs, with respect to the related issue of fragmentation. 
Here we follow up the work of Rice et al., in order to quantify under which conditions a 
local model for transport in self-gravitating discs is justified. 



NUMERICAL SETUP 

The simulations presented here were performed u sing smoothed particle hydrodynam- 
ics (SPH), a Lagrangian hydrodynamics code [see[lj]. Our code uses a tree structure to 
calculate gravitational forces and the nearest neighbours of particles. The particles are 
advanced with individual time-steps [ 14], resulting in an enormous saving in computa- 
tional times when a large range of dynamical timescales are involved. 

We consider a system comprising a central star, modelled as a point mass M^, sur- 
rounded by a disc with mass M^u^c. We have performed several simulation using different 
values of the parameter q = Mdisc/Mi,. The initial surface density profile was taken to be 



Figure 1. Equatorial density structure at the end of the simulations for (left) q = 0.05, (middle) ^ = 0.1, 
and (right) ^ = 0.25. 



a power-law E . The disc is initially gravitationally stable and the minimum value 
of Q is Qmin = 2. Our calculations are essentially scale-free. In dimensionless units, the 
disc extends from = 0.25 to Vout = 25. We have used = 250000 particles. 

We use an adiabatic equation of state, with adiabatic index 7= 5/3. We explicitly 
solve the energy equation for our system, including heating from shocks and from 
artificial viscosity and introducing a cooling term, simply parameterized as: 

f^) =--i^L. (3) 

V / cool ^cool 

where m, is the internal energy of a representative particle, and, as in Gammie lHoll and 
in Rice et al. [12|], the cooling time-scale tcooi = K where /3 is a free parameter, 
independent of radius. Gammie [10] and Rice et al. flVM have shown that if cooling 
is too fast, then the disc will fragment as a result of the gravitational instabilities. In 
our simulations we have used /3 = 7.5, in which case none of our discs was found to 
fragment. 



RESULTS 

We have run our simulations for a few thermal timescales, in order to reach thermal 
equilibrium. In all cases the disc initially cools down until the value of Q becomes small 
enough for gravitational instabilities to become effective and to provide a source of 
heating to the disc. At later stages the disc settles into a self-regulated state with 2 ~ 1 
over a large portion of the disc. 

The component of the stress tensor (integrated in the z-direction) associated with 
self-gravity is given by 

rj^grav _ f J Mi (4\ 




Figure 2. Effective a produced by gravitational instabilities for (left) q = 0.05, (middle) ^ = 0.1, and 
(right) q — 0.25. The top panel shows the separate contribution of Ugrav and UReyn, the lower panel shows 
the sum of the two contributions compared with the expected value from a local viscous model (dotted 
line). 



where g is the gravitational field of the disc. The hydrodynamic (Reynolds) contribution 
to the stress tensor is given by: 



r,Reyn 
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where 5v = v — u is the velocity fluctuation, v is the fluid velocity and u is the mean 
fluid velocity. 

After averaging the stress tensor azimuthally and radially, over a small region Ar = 
O.lr, we have computed the corresponding value of a (see Eq. ([T])): 
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The resulting radial profiles of a are shown in Fig.|2lfor the three cases q = 0.05, q = 
0. 1 , and ^ = 0.25. The upper panels show separately the hydrodynamic and gravitational 
contributions to a, while the bottom panels show the sum of the two. The results are also 
averaged in time, over the last 500 timesteps of the simulations. 

The simple parameterization of the cooling time adopted here allows us to have a 
direct insight on the problem of characterizing transport properties and dissipation in the 
disc, by using the following result [16], which follows directly from the energy balance 
equation: if the disc is in thermal equilibrium, and if heating is dominated by a purely 
viscous process, described by Eq. dU), then the viscosity coefficient a and the cooling 
time-scale t^ooi satisfy the following relation: 
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Figure 3. Contribution to the stress at r = 15 from regions of the discs at a distance d from r. Dashed 
hne: q — 0.05, solid line: ^ = 0.1. Dotted line: q = 0.25 

In our simulations tcooi = T-5Q.^\ so from Eq. Q we would expect that in thermal 
equilibrium, if energy dissipation can be described locally, a should be approximately 
constant with radius (this is indicated with the dotted line in the bottom panels of 
Fig. 121). We can interpret the bottom panels of Fig. |2l in the following way: the solid 
line gives a measure of the actual torque induced by gravitational instabilities in our 
simulations, while the dotted line is the torque that a viscous process should exert in 
order to dissipate the amount of energy which is actually dissipated in the simulations. 
Therefore, when the computed value of a is larger than the expected one, gravitational 
instabilities dissipate less energy than a viscous process, whereas in the opposite case, 
they dissipate more energy than a viscous process. Our results indicate that, up to a disc 
mass Mjisc = 0.25M^, gravitationally induced transport is fairly well described within a 
local framework. This result could also be anticipated, based on Fig.[Tl that shows how 
the disc dynamics is dominated by high-m modes, that dissipate on a short length-scale. 

As a separate test for the locality of transport, we have also computed apart {r,d), 
which is defined as the gravitational part of a(r), where the gravitational field g is 
computed taking into account only those particles which are inside a spherical radius 
d from the radial point r at which we are computing the stress. This quantity gives us a 
measure of the size of the region that mostly contributes to the stress at a given point. 

Fig. m shows the results for r = 15 for the three simulations. The upper panel shows 
clearly that the more massive the disc is, the more far regions contribute to the stress. 
However, it should be noted that, since 2 ~ 1, a more massive disc is also hotter, and 
therefore thicker. The bottom panel of Fig. El shows OCpart as a function of d/H, where H 
is the disc thickness. In all cases, more than 80% of the contribution to the gravitational 
torque come from a region with size A ^ lOH. Transport is going to be local if A ^ r. 
We can therefore conclude that, as long as Mjisc < 0.25A/*, a disc with H/r <ti 0.1 will 
be reasonably described in terms of a local model. This is basically the global analogous 
of the results obtained by Gammie lHoll based on a local, 2D model. 




Note that detailed vertical structure models of FU Ononis objects |5] lead to disc 
thickness H /r ^O.l, so that in this case we could expect global effects to be important 
in determining the emission properties of the disc, as suggested by Lodato and Bertin 



CONCLUSIONS 

We have performed global, three-dimensional SPH simulations of self-gravitating ac- 
cretion discs. Our results show how the heating provided by gravitational instabilities 
can balance the cooling (that we have parameterized in a simplified way), leading to a 
self-regulated state where 2 ~ 1. With a cooling rate tcooi = 7.5t2^\ none of our sim- 
ulations was found to fragment. We have characterized the transport properties of the 
disc, by computing the stress tensor associated with gravitational instabilities and com- 
paring the corresponding viscous dissipation rate to the actual heating rate in the disc. 
We have found that, if the disc is less massive than 0.25M*, gravitational disturbances 
are dominated by high-m modes, so that the transport properties of the disc can be de- 
scribed reasonably well in terms of a local, viscous formalism. In particular, we have 
found that the gravitational torque comes from a region of size ^ \0H, where H is the 
disc thickness. We can therefore expect that discs with an aspect ratio H/r <0.l can be 
treated locally. 
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